/***************************************************************************//**
 * @file MSc_IntegratorRungeKutta4.h
 * @author  Marek M. Cel <marekcel@mscsim.org>
 *
 * @section LICENSE
 *
 * Copyright (C) 2013 Marek M. Cel
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 ******************************************************************************/
#ifndef MSC_INTEGRATOR_RUNGE_KUTTA_4_H
#define MSC_INTEGRATOR_RUNGE_KUTTA_4_H

////////////////////////////////////////////////////////////////////////////////

#include "MSc_Integrator.h"

////////////////////////////////////////////////////////////////////////////////

/**
 * Runge-Kutta 4th order numerical integration template class.
 * @see Krupowicz A.: Metody numeryczne zagadnien poczatkowych rownan rozniczkowych zwyczajnych (polish)
 * @see Baron B., Piatek L.: Metody numeryczne w C++Builder (polish)
 */
template <class TYPE>
class MSc_IntegratorRungeKutta4 : public MSc_Integrator<TYPE>
{
public:

    /** Constructor. */
    MSc_IntegratorRungeKutta4( int size = 0, TYPE *obj = 0, void (TYPE::*fun)(const double *, double *) = 0 ) :
        MSc_Integrator<TYPE> ( size, obj, fun ),

        m_k_1 ( 0 ),
        m_k_2 ( 0 ),
        m_k_3 ( 0 ),
        m_k_4 ( 0 ),

        m_x_t ( 0 )
    {
        m_k_1 = new double [ this->m_size ];
        m_k_2 = new double [ this->m_size ];
        m_k_3 = new double [ this->m_size ];
        m_k_4 = new double [ this->m_size ];

        m_x_t = new double [ this->m_size ];
    }

    /** Destructor. */
    virtual ~MSc_IntegratorRungeKutta4()
    {
        if ( m_k_1 ) delete [] m_k_1; m_k_1 = 0;
        if ( m_k_2 ) delete [] m_k_2; m_k_2 = 0;
        if ( m_k_3 ) delete [] m_k_3; m_k_3 = 0;
        if ( m_k_4 ) delete [] m_k_4; m_k_4 = 0;

        if ( m_x_t ) delete [] m_x_t; m_x_t = 0;
    }

    /**
     * Integrates given vector using
     * Runge-Kutta 4th order integration algorithm.
     * @param step integration time step [s]
     * @param vect integrating vector
     */
    void integrate( double step, double *vect )
    {
        for ( int i = 0; i < this->m_size; i++ )
        {
            m_x_t[ i ] = vect[ i ];
            m_k_1[ i ] = m_k_2[ i ] = m_k_3[ i ] = m_k_4[ i ] = 0.0;
        }

        // k1 - derivatives calculation
        this->fun( m_x_t, m_k_1 );

        // k2 - derivatives calculation
        for ( int i = 0; i < this->m_size; i++ ) m_x_t[ i ] = vect[ i ] + step * m_k_1[ i ] / 2.0;
        this->fun( m_x_t, m_k_2 );

        // k3 - derivatives calculation
        for ( int i = 0; i < this->m_size; i++ ) m_x_t[ i ] = vect[ i ] + step * m_k_2[ i ] / 2.0;
        this->fun( m_x_t, m_k_3 );

        // k4 - derivatives calculation
        for ( int i = 0; i < this->m_size; i++ ) m_x_t[ i ] = vect[ i ] + step * m_k_3[ i ];
        this->fun( m_x_t, m_k_4 );

        // integration
        for ( int i = 0; i < this->m_size; i++ )
        {
            vect[ i ] = vect[ i ]
                      + step * ( m_k_1[ i ] + 2.0 * m_k_2[ i ] + 2.0 * m_k_3[ i ] + m_k_4[ i ] ) / 6.0;
        }
    }

private:

    double *m_k_1;  ///< auxiliary vector
    double *m_k_2;  ///< auxiliary vector
    double *m_k_3;  ///< auxiliary vector
    double *m_k_4;  ///< auxiliary vector

    double *m_x_t;  ///< auxiliary vector

    /** Using this constructor is forbidden. */
    MSc_IntegratorRungeKutta4( const MSc_IntegratorRungeKutta4 & ) {}
};

////////////////////////////////////////////////////////////////////////////////

#endif // MSC_INTEGRATOR_RUNGE_KUTTA_4_H


